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We present momentum resolved single-particle spectra for the three-dimensional Hubbard model 
for the paramagnetic and antiferromagnetically ordered phase obtained within the dynamical cluster 
approximation. The effective cluster problem is solved by continuous-time Quantum Monte Carlo 
simulations. The absence of a time discretization error and the ability to perform Monte Carlo 
measurements directly in Matsubara frequencies enable us to analytically continue the self-energies 
by maximum entropy, which is essential to obtaining momentum resolved spectral functions for the 
Neel state. We investigate the dependence on temperature and interaction strength and the effect 
of magnetic frustration introduced by a next-nearest-neighbor hopping. One particular question we 
address here is the influence of the frustrating interaction on the metal-insulator transition of the 
three-dimensional Hubbard model. 

PACS numbers: 71.10.Fd 71.15.-m 71.28.+d 71.30.-h 



I. INTRODUCTION 

One of the paradigms for correlation effects and com- 
peting orders in solid state physics is the Hubbard 
modePl2l 




The operators [cia) create (annihilate) an electron 
with spin cr £ {tii} at lattice site i, riia = cl^Cia is 
the particle number operator, t describes the hopping 
between neighboring sites (denoted by (i, j)), t' the hop- 
ping between next-nearest neighbors (denoted by ((z, j))) 
and U implements the local Coulomb repulsion. 

The Hubbard model, despite its simple structure, can 
only be solved exactly in on^ and infinit^ spatial di- 
mensions. There are several expectations one can de- 
duce from general energetic arguments and in partic- 
ular from the connection between the Hubbard model 
and the Heisenberg hamiltonian in the limit of large in- 
teraction strength t/U — > 0.^ For a three-dimensional 
simple-cubic lattice and t' — 0, the Hubbard model at 
half filling shows antiferromagnetic order at finite tem- 
perature for any value U > 0. The doubling of the unit 
cell causes this ordered state to be an insulator. With 
increased next-nearest-neighbor hopping this transition 
is suppressed and a Mott-Hubbard metal-insulator tran- 
sition (MH-MIT) is expected to appear in the param- 
agnetic state at some non-zero critical value Uc of the 
interaction. 

Quantum Monte Carlo (QMC) methods are powerful 
tools that enable the controlled calculation of proper- 
ties of large interacting quantum many-particle systems. 
Examples include spin models and many bosonic sys- 



tems. However, simulations of fermionic models away 
from particle-hole symmetry are often severely hampered 
by the fermionic sign problem.^ In particular, the identi- 
fication of ordered phases, which requires a reliable finite- 
size scaling, becomes exceedingly complicated. Indepen- 
dent of the sign problem, the direct investigation of the 
properties of ordered phases possibly present in the ther- 
modynamic limit is not possible, since any QMC simula- 
tion is performed on a finite system which cannot exhibit 
a spontaneously broken symmetry. 

Therefore an approximation scheme allowing (i) cal- 
culations in the thermodynamic limit while (ii) includ- 
ing dynamical correlations in a controlled way is highly 
desirable. The dynamical mean-field theory (DMFT)^'^ 
and its cluster extensions^ are such theories. The DMFT 
maps the lattice problem onto an effective single-site im- 
purity model, at the cost of neglecting non-local many- 
body correlation effects. In many cases, the sign prob- 
lem of the resulting impurity model is either absent or 
manageable. However, these non-local correlation effects 
are often crucial for the interplay between Fermi liquid 
and more exotic states of matter. Cluster mean-field 
theories are extensions of the DMFT to finite clusters, 
re-introducing non-local (short-ranged) correlations in a 
systematic manner but at the same time increasing the 
complexity of a simulation. In the limit of infinite cluster 
size they become exact. The controlled extrapolation of 
cluster results to th e thermodynamic limit is often feasi- 
ble in practice! ^° * ^^ * 

DMFT has proven to be a very powerful tool for study- 
ing the fundamental aspects of the MH- MF#I. Within 
the DMFT, however, the MH-MIT is completely hidden 
inside an antiferromagnetic phase, which is insulating by 
symmetry.'^ Introducing magnetic frustration, e.g., by 
a non-zero t', tends to suppress the magnetic order and 
shift the critical interaction strength Uc toward lower val- 
ues. Thus the MH-MIT eventually emerges from the an- 
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tiferromagnetic phase for large values of 

In low-dimensional systems (single site) DMFT is in 
general not a good approximation. In particular, for the 
one-dimensional Hubbard model, non-local correlations 
are in fact dominant)^ leading to a complete breakdown 
of Fermi liquid physics and the formation of a novel low- 
energy fixed point, the Luttinger liquidP Similarly, in 
two dimensions a strong influence of spin fluctuations 
in the Hubbard model is expected, in particular at and 
close to half filling. Since the Mermin- Wagner theorem 
forbids the formation of an ordered state in two dimen- 
sions at finite temperature, the existing strong magnetic 
correlations will lead to correspondingly strong dynam- 
ical fluctuations at T > and will possibly trigger a 
similar breakdown of Fermi liquid physics as in the one- 
dimensional case. Evidence for this behavior h as in deed 
been observed in various numerical simulations.'^'^^'^^ 

The importance of short-ranged correlations in the 
three-dimensional Hubbard model is less clear and less 
well studied, and detailed studies of the phase dia- 
gram at hig h temperature have only recently begun to 
appear.fiSESl On the one hand, the precise value of the 
Neel temperature and critical exponents for the tran- 
sition into the antiferromagnetic state will be directly 
influenced by the presence of spin fluctuations.'^ Since 
the antiferromagnetically ordered phase is an insulator, 
we expect the antiferromagnetic spin fluctuations in the 
paramagnetic phase to stabilize the MH-MIT, thus shift- 
ing the critical Uc toward lower values. Adding frustra- 
tion by, e.g., next-nearest-neighbor hopping t' will fur- 
ther enhance this effect. One may surmise that the MH- 
MIT will eventually emerge from the antiferromagnetic 
phase as in DMFT. To our knowledge this has not yet 
been investigated in detail. 

In this work we study the Hubbard model for a three- 
dimensional cubic l attice using the dynamical cluster ap- 
proximation (DCA)28I291 cluster dynamical mean-fleld al- 
gorithm on clusters of size 18. We present momentum 
resolved single-particle spectra in the paramagnetic and 
in the antiferromagnetic phase, and investigate the in- 
fluence of frustration effects caused by a next-nearest- 
neighbor hopping t' . We focus our investigation on the 
interplay of frustration and spin fluctuations in the vicin- 
ity of the paramagnetic metal insulator transition. 



II. METHOD 

We study the Hubbard model in three dimensions 
within the DCA to include both the short- to medium- 
ranged antiferromagnetic fluctuations and the possibility 
of actual long-range antiferromagnetic order. Since the 
DCA maps the lattice problem onto an effective peri- 
odic cluster coupled to a dynamic bath, numerically ex- 
act quantum Monte-Carlo (QMC) algorithms are ideally 
suited to solving this effective model. 

Of particular interest in correlated electron systems are 
dynamical correlation functions such as single-particle 



spectra. However, QMC provides data only on the imag- 
inary time or frequency axis, and the necessary analytic 
continuation of these data has proven to be difScult. The 
standard tool used to solve this problem is the maximum 
entropy method (MEM).3o 

Previously, the quasi-standard for simulations of 
fermionic many-particle systems was the Hirsch-Fye 
algorithm, ■^^ which uses a discretization of the imaginary 
time axis. An alternative has evolved in recent years 
by the development of QMC algorithms in continuous 
imaginary time.[2^^ The absence of a time discretiza- 
tion error and the possibility of Monte Carlo measure- 
ments directly in Matsubara frequencies^ enhance the 
quality of the data significantlj^^ and hence enable us to 
directly analytically continue self-energies.'^ This avoids 
the extraction of the self-energies from already contin- 
ued Green's functions by a numerically difficult multi- 
dimensional root finding algorithm .I^Sl In this paper we 
use an implementation of the continuous-time interaction 
expansion (CT-INT) QM C algo rithm initially described 
by Rubtsov and co-workers^^ES which performs a system- 
atic expansion in the interaction term of the Hamiltonian. 

For a simple cubic lattice in three dimensions the 
dispersion including nearest- and next-nearest- neighbor 
hopping reads 



efe ~ — 2t cos(fci) — At' [cos(fci) cos(fc2) 



-C0s(fc2) C0s(fc3) -I- COs(fci) COs(fc3)] , (2) 



where k = (ki, fe, ^3) is an element of the first Brillouin 
zone of the simple cubic lattice. The full bandwidth 
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12t for \f \ < t/A 

8t + 16|i'| iov\t'\>t/A 



(3) 



of this dispersion is used as the energy scale in this paper. 

Let us briefly recall the essential aspects of the DCA. 
The central quantity is the single particle Green's func- 
tion in imaginary time t, defined by 

G..,(r) = -(rc.,(r)4.) . (4) 

Here T is the imaginary-time ordering operator, (•) 
denotes a thermal expectation value, and Caiij) — 
The spatial and temporal Fourier trans- 
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form of Green's function is 



X y dr exp (ia;„r) G„ij{j) 



(5) 



where k is located in the first Brillouin zone and w„ = 
(2n -I- l)7r//3 with n G Z and /3 Xjk^T denotes the 
fermionic Matsubara frequencies. T is the temperature. 
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FIG. 1: (Color online) Dynamical cluster approximation il- 
lustrated for an eight-site cluster in two dimensions. In real 
space (a) the origin of a cluster is labeled by r. Each site of 
the cluster is identified by R. A Fourier transformation maps 
the coordinate r = r + R oi each lattice site to a vector fc in 
the first Brillouin zone (b). The cluster momentum K now 
identifies the center of a cell in momentum space. All points 
inside this patch can be reached by k. The DCA integrates 
out k and thus replaces the full fc dependence of the lattice 
by the cells labeled by K. 



and fee the Boltzmann's constant. With Dyson's equa- 
tion, we write Green's function as 
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(6) 



thus introducing the single-particle self-energy 'E^ki'i-^n) 
which contains all the many-body correlation effects in 
the system and, in general, is a function of both momen- 
tum k and energy w„. 

Within DCA, the full lattice model is approximated 
by a finite cluster of size N embedded in a mean field. 
We tile the first Brillouin zone into N non-overlapping 
cells, each represented by its central momentum K [see 
Fig. [ijb) for an example] . The full fe-dependence of the 
model is approximated by the discrete set of N cluster 
momenta K by setting Eo-fc(iw„) « Eo-i<:(ia;„). Averag- 
ing over the volume V of the cell corresponding to cluster 
momentum K, one obtains the quantity 



dk, 
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, (7) 



which defines an effective non-interacting cluster via 

^?<Tjt(iw„)"^ = Gcta: + So-/<:(iw„) . (8) 



QMC Cluster Solver 
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Gk = y I Ak [lujn + /i - e^+fe 



FIG. 2: Self-energy Eo-K:(iti;„) is determined self-consistently 
by iterating the depicted procedure until convergence is 
reached. The bottom line shows the calculation of the coarse- 
grained Green's function Gk by averaging over the momen- 
tum patch centered around K via integrating fc over the vol- 
ume V of the patch. The dependency of Green's functions 
and the self-energy on iujn and a is omitted for simplicity. 



With this set of quantities, a suitable method to solve 
the effective cluster defined by QaKi^'^n) and the inter- 
action J7, one can determine the self-energy ScrK:(iw„) 
iteratively as depicted in Fig. [2j 

From the QMC algorithm used to solve the effective 
cluster, we obtain the cluster Green's function Go-K"(iwn}. 
Usually one then uses the maximum entropy methocP^l 
to analytically continue this quantity to the real axis. 
In order to be able to reverse the coarse-graining, i.e., 
calculate Gcrk,{<^ + iO^) for all fc from the first Brillouin 
zone, one needs access to the self-energy So- a: 
which needs to be obtained by numerical inversion of 
Eq. [7] While this is feasible in the paramagnetic phase, 
the matrix structure appearing in the antiferromagneti- 
cally ordered phase (see section IV I renders this approach 
impractical. 



We follow here an alternative route and analytically 
continue the self-energy insteacPH, which is related to the 
cluster Green's function by 



(9) 



i.e., an inversion of Ga-Ki^i-^n)- This inverse is calcu- 
lated directly from the Monte Carlo bins using a jack- 
knife procedur^^ and therefore incorporates a full error 
propagation of the covariance matrix. The bare Green's 
function Ga-Ki^'-^n) is viewed here as an input parameter; 
and error propagation of errors contained in our estimate 
of it, which would require error propagation over sub- 
sequent iterations, is not considered here, such that all 
errors in GaKi^^n) are neglected. 

The analytic continuation of the self-energy from imag- 
inary to real frequencies is then performed by the max- 
imum entropy method,'22! using a standard implementa- 
tion of the algorithm following^. To accurately continue 
self-energies with MEM, their high frequency behavior 
has to be known"^^. To this end we perform a high- 
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frequency expansion of the self-energy 

^.Ki^^n) = + ^ + , (10) 

where the coefficients are given by (see Appendijj]) 

K^U [{n.„) ^) (11) 



and 

El =[/2(n_.)(l-(n_.)) 
We now define the quantity 



S^x(iw„) := 



E^i<:(icJn) - £° 
Si 



(12) 



(13) 



Since the average number density (n_o.) is a Monte Carlo 
measurement, we estimate S^^^ (iaj„) and its covariance 
matrix by a jackknife procedure. The rescaled self-energy 
Ti'^j^iiUn) as function of Matsubara frequencies is related 
to the imaginary part lTiiY,'^j^{uj + iO^) on the real fre- 
quency axis through the Hilbert transform 



dc.'^^^ 



{J) 



(14) 



By virtue of the rescaling Eq. [13] we furthermore have 

oo 

j dujlmY.'„K{^)^\, (15) 

— oo 

i. e., the spectral function — IrlmE'^j^ (w) is non-negative, 
normalized to unity, and can thus be calculated by the 
MEM from the data on the imaginary axis. The real part 
of the self-energy then follows from the Kramers-Kronig 
relation 



oo 

ReE.^(.) = -ip/dc.'^^^^^^ + EO, 



(16) 



where P J denotes a principal value integral. An example 
for a full self-energy on the real-frequency axis is shown 
in Fig. [Sf^a). An interpolation of the coarse-grained self- 
energies yields the self-energy Yi„]^{uj) for all momenta 
k of the Brillouin zone. We use a three-dimensional in- 
terpolation based on Akima splineJ^ which provide a 
smooth interpolation along the momentum points while 
avoiding spurious oscillations. Finally, the single-particle 
spectral function A^ki'jj) is calculated using Dyson's 
equation: 



— Ini — 



- /i - efc - Scrfc(w) 



(17) 



The quantity A„f.(io)co\i\A in principle also be calcu- 
lated by an analytic continuation of the Green's function 
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FIG. 3: (Color online) (a) Self-energy on the real-frequency 
axis in the paramagnetic phase for K = T , U = W, t' = 0, 
and r = 0.02 (a). The default model that entered the 
MEM calculation of the imaginary part is also shown. The 
real part is obtained from the imaginary part via Eq. |16| 
(b) Real and imaginary parts of the interpolated self-energy 
Sfe(ia;o) for the smallest Matsubara frequency for an 18-site 
cluster, U = 0.67 W and T = 0.03 W at half filling. The 
self-energy values at momenta between the discrete cluster 
momenta are obtained from an interpolation using Akima 
spline in three dimensions. The horizontal straight lines 
denote DMFT results. The interpolation follows the path 
along the high symmetry points of the first Brillouin zone de- 
picted in panel (c). For an estimate of finite-size effects in (b) 
see also Ref.'^ 
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FIG. 4: (Color online) Momentum resolved single-particle spectra for T = 0.02 W and t' — 0. The momenta k follow a path 
along high symmetry directions of the first Brillouin zone as depicted in Fig. |3|;. The left part of each diagram shows the local 
single-particle spectrum A{ijj) derived from the direct analytic continuation of the Green's function. The dashed line denotes 
the bare dispersion efc. 



at each cluster momentum K followed by an interpola- 
tion of the Green's function values directly. However, this 
procedure approximates the exactly known momentum 
dependence of the dispersion e^, which causes a signifi- 
cant loss of momentum resolution. We will not consider 
it any further. 

As an example, Fig. [sjjb) shows the interpolated 
momentum-resolved self-energy for the smallest Matsub- 
ara frequency iloq (full and dashed lines). The results 
were obtained for a cluster with 18 K momenta and a 
Coulomb repulsion U = 0.67 W at T = OMW. The 
comparison with the corresponding DMFT result (dotted 
lines) shows that the momentum resolution of the DCA 
adds significant k dependence to the self-energy. Thus 
the many-particle renormalizations acquire a significant 
fc dependence in DCA, even for the three-dim ensio nal 
Hubbard model at comparatively weak coupling .1^2111] 

The available computational resources and the qual- 
ity of data needed for high-precision analytic continua- 
tion, limit us to study clusters of comparatively small 
size. We limit ourselves to a study of a cluster of size 
= 18 described by the basis vectors ai = (1,1,2), 
a2 = (2, 2, -2), and = (2, -1, -l)^his cluster is the 
optimal bipartite cluster of this size^ following the cri- 
teria proposed by Belts et al^^. Since we are primarily 
interested in identifying trends and basic physical effects, 
we do not perform calculations for larger clusters to ob- 
tain a finite-size scaling as would have been necessary. 



e.g., for a precise estimation of the equation of st ate o r 
the Neel temperature in the thermodynamic limit .1^21211 

III. PROPERTIES OF THE PARAMAGNETIC 
PHASE 

We begin the discussion of our results by present- 
ing spectral functions in the paramagnetic phase of the 
model, i.e., we manually suppress long-range order. This 
allows us to distinguish dynamical effects coming from 
fluctuations from effects caused by the (static) symme- 
try breaking. 

A. Metallic phase 

Figure |4] presents single-particle dispersions in the 
paramagnetic phase for T = 0.02 W, t' = and four dif- 
ferent Coulomb repulsions U = 0.67W, U = 0.8W, U = 
0.87W, and U = W. The selected k points follow a 
path along the high symmetry points of the first Bril- 
louin zone depicted in Fig. [sjjc). Local single particle 
spectra are shown on the left of each panel. For small 
U = 0.67 W we observe a clear quasi-particle peak at the 
Fermi level, both in the density of states (DOS) and the 
spectral function. The momentum resolved spectra show 
that the main contributions to the quasi-particle peak 
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FIG. 5: (Color online) Momentum resolved single-particle 
spectra for T = 0.02 W, t' = 0, and U = 0.93 W using a 
non-dispersive DMFT self-energy (a) and a momentum de- 
pendent cluster self-energy (b). The inset highlights a part of 
the spectrum using an alternative color scheme. 



are situated halfway between the T and R points and the 
X and M points, respectively. Comparison of the peak 
dispersion in these regions to the non-interacting disper- 
sion (dashed line) reveals a clear flattening at the Fermi 
energy, i.e., an increased effective mass of the quasi- 
particles. At higher energies additional structures — the 
lower and upper Hubbard bands — are visible. They fol- 
low the curvature of the bare dispersion but are shifted 
to higher energies and connected to the quasi-particle 
band through broad "waterfall" -like pieces reminiscent 
of structures observed in angle-resolved photoemission 
spectroscopy of cuprates.^ For increasing Coulomb re- 
pulsion U the quasi-particle band at the Fermi energy 
vanishes and is replaced by an insulating gap. At the 
same time the dispersion of the high-energy structures 
flattens: the system becomes more localized. Thus for 
the temperature studied a crossover from a metallic dis- 
persion ai U = 0.67 W to a Mott insulating phase at 
U = W is clearly visible. We will return to the details of 
the Mott-Hubbard metal-insulator transition in section 

cnBi 

To make the influence of the momentum dependence 
of the self-energy in the spectra more transparent. Fig. [5] 



compares an insulating spectrum based on self-energies 
of the 18 site cluster to the corresponding spectrum based 
on a momentum-independent DMFT self-energy. While 
many overall features are similar, there are qualitative 
differences. For example, a blowup of the details of the 
spectra close to the F point (insets to Fig. [s]) shows that 
a substantial part of the DMFT spectrum around the F 
point is located just above the Fermi energy. This con- 
tribution is shifted to higher frequencies in the cluster 
calculation, and the curvature is reversed, more resem- 
bling the structure of the lower Hubbard band but with 
much less spectral weight. The feature can be interpreted 
as a precursor of the complete symmetry with respect 
to the Fermi energy occurring for spectra in the anti- 
ferromagnetically ordered phase (see section IV). Thus, 



we attribute these pale reflections of the lower Hubbard 
band to the so-called shadow hand^^ arising due to an- 
tiferromagnetic fluctuations not contained in the DMFT 
simulation. 



Next we examine the influence of a next-nearest neigh- 
bor hopping t' . In quantum Monte Carlo calculations, 
a non-zero value of t' introduces a fermionic sign prob- 
lem even at half filling that may lead to a significantly 
larger computational cost. However, for the tempera- 
tures, Coulomb repulsions and cluster size studied here 
the average sign is always greater than 0.94 and thus 
affects the efficiency of the simulations only marginally. 

Fig.[6]shows results of calculations for U = 0.67 W and 
different values for t' . The particle-hole symmetric spec- 
trum of Fig. |4][a) becomes more and more asymmetric 
with increasing t' . We attribute most of these changes 
to the changes in the bare dispersion e^. However, while 
for small to moderate t' the quasi-particle properties do 
not seem to change much, one observes a significant re- 
duction in the spectral weight at the Fermi energy for 
larger t' , resulting in a reduction of the quasiparticle 
peak in the analytically continued spectra. For exam- 
ple, at t' = —t the integrated weight of the quasi-particle 
peak is reduced by 50% compared to the value for t' ^ 
in Fig. Qa) . These observations indicate a reduction of 
the quasi-particle mass with increasing t' , in accordance 
with previous DMFT finding^. 

Frustration effects on the insulating spectrum in 
Fig. ^d) {U = W) are illustrated in Fig. [t] As for the 
metallic spectrum, the features present in Fig. |4]Jd) for 
t' = initially change only weakly and in particular the 
"shadow" structures seem to be present, albeit with re- 
duced weight. For strong frustration, Hubbard bands 
become dominant. These Hubbard bands have a rather 
well-defined structure and dispersion reminiscent of the 
bare dispersion. Furthermore, a peak develops just below 
the Fermi energy, which becomes more pronounced with 
increasing t' and appears to have only small momentum 
dependence. 
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FIG. 6: (Color online) Single-particle spectra for U = 0.67 W, T = 0.02 and non-zero t' . 
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FIG. 7: (Color online) Single-particle spectra for U = W, T = 0.02 W and non-zero t' . 
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B. Mott-Hubbard Metal-Insulator Transition 

One of the interesting properties of the Hubbard model 
is the formation of a correlation driven metal-insulator 
transition in the paramagnetic phase, the so-called Mott- 
Hubbard metal-insulator transition (MH-MIT) . Different 
from the conventional band or Slater insulators for even 
electron number, where the insulating behavior is due to 
a completely filled band, the MH-MIT occurs in a par- 
tially filled band, which within a simple single-particle 
picture would thus be conducting. Such a transition 
is believed to frequently be present in transition metal 
oxides,'^ and has been discussed both from an experi- 
mental and theoretical point of view over the past 50 
years. Note that it is also different from the insulating 
state originating from an antiferromagnetic order; the 
latter implies a broken translational invariance and can 
hence be purely explained by band structure effects (see 
also section IV). 



On a simple cubic lattice with only nearest-neighbor 
hopping the MH-MIT of the Hubbard model at half 
filling is completely covered by the antiferromagnetic 
phase.!^ Nevertheless, one may study it within a general- 
ized mean-field theory by suppressing long-range antifer- 
romagnetic order in the system. This is accomplished 
by enforcing translational symmetry on the bath and 
thereby preventing any symmetry-breaking. One of the 
clearest signs, numerically, for identifying the MH-MIT 
is given by the value of the DOS at the Fermi energy: 
A non-zero value at T = is indicative of a metal, a 
zero value of an insulator. Identifying the MIT at T > 
without a detailed analysis of the temperature depen- 
dence is more subtle, but again the DOS at the Fermi 
level can serve as "order parameter" : Away from the 
MH-MIT the DOS varies smoothly as function of tem- 
perature. When approaching the MH-MIT at half filling 
as function of U or T, DM FT analyses suggest that the 
DOS will show a jumpSHSbelow a critical end point. Fur- 
thermore, the transition is of first order with a nice hys- 
teresis in the critical region. This behavior has been con- 
firmed for two-dimensional systems within cluster DMFT 
on smalP^and largei'^^ clusters. Note that we do not dis- 
cuss the approach to the MH-MIT as function of doping. 
The behavior of the system at this transition is very dif- 
ferent, and has been addressed by a number of groups 

both in DMFT and cluster variants on a wide range of 
systems.lSlElHIIl 

We estimate the density of states at the Fermi 
level from the low-frequency behavior of our Matsub- 
ara Green's function Gii{iuJn) that is available as direct 
Monte Carlo measurement and does not require analytic 
continuation. This approach is based on the relation 



cosh ^uj 



-ttAkIO) (18) 



for the long-time behavior of the imaginary time Green's 
function and the fact that the long-time behavior trans- 
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FIG. 8: (Color online) Imaginary part of the local Green's 
function in Matsubara frequencies for t' = 0, T = 0.01 W 
(a) and T — 0.0151^ (b). Several values of U around the 
MH-MIT are shown. 



lates into the low-energy behavior under the Fourier 
transform. 

Figure [8](a) shows the imaginary part of Gioci^^n) '■— 
Gii{iu}n) for t' — and T = 0.01 W obtained within the 
DCA for a cluster size of 18. The jump from an insulating 
Green's function at U = 0.80 to a metallic solution at 
U = 0.766 clearly shows the location of the MH-MIT 
at this temperature. For T = 0.015 we could only de- 
tect a crossover from insulating to metallic behavior [see 
Fig.jsj^b)]. This indicates that the critical endpoint of the 
MH-MIT transition hue is located between T = 0.01 W 
and T ~ 0.015 VF, substantially below the Neel temper- 
ature at t his in teraction strength [Tn = 0.030(3) W at 
U = 0.8^^^. 

Another observable that shows a clear signal of the 
MH-MIT is the effective mass, defined as 



'"-fc 
m 



= 1 



(19) 



where m denotes the bare carrier mass. The effective 
mass at non-zero (but sufficiently low) temperature may 
be estimated from the QMC without resorting to analyt- 



9 



t' = -0.2 1, U = 0.73 W 
t' = -0.2 1, U 0.80 
<' = 0, [/ = 0.73 W 
i' = 0, [/ 0.80 




FIG. 9: (Color online) (a) Quasi-particle mass as func- 
tion of mean distance d from the Fermi surface in units of 
the lattice spacing a (a). The distance djc is defined as the 
average of the distances between the Fermi surface and all mo- 
menta in the cluster cell described by the cluster momentum 
K. The 18 cluster momenta reduce to five different masses 
due to point symmetries. The labels 1 to 5 refer to the points 
depicted in panel (b), which also shows the non- interacting 
Fermi surface for t' = 0. 



ical continuations 



m 



ImEcrfc(iwo) 



UJQ 



(20) 



where ujq is the lowest Matsubara frequency (see, e.g., 
Fig. 10 in Ref.l^ for a comparison between directly ob- 
tained and analytically continued estimates). At the MH- 
MIT, across the Fermi surface exhibits a sharp in- 
crease. 

The estimate for the effective mass obtained that way 
is shown in Fig. ^a) for two different values U = 
0.73W < Uc{t' = 0) and [/ = 0.8 > Uc{t' = 0) and 
t' = and -0.2 i at T = 0.01 VF. We do not show the 
result for from the interpolated data, as the division 
by Wo with oJo/W <^ 1 also strongly enhances spurious 
artificial oscillations due to interpolation, but rather we 
present the masses m*j^ for each of the 18 cluster mo- 
menta as functions of their mean distance to the non- 
interacting Fermi surface for t' = 0. The mean distance 



is calculated by averaging the distance to the Fermi sur- 
face of all points inside the cluster cell around K. Due 
to symmetry, some cluster momenta are equivalent and 
thus we obtain only five different effective masses. Fig- 
ure |9]^b) depicts one representative cluster momentum for 
each one of these five equivalence classes. 

For both values of t' the K dependence of m*j^ for U 
deep inside the metallic phase is rather weak, although 
nevertheless visible. The insulating phase, in contrast, 
shows a dramatic increase of the effective mass for the 
cluster momenta near the Fermi surface. For example, 
point 1 with the strongest enhancement is situated on the 
Fermi surface of the non-interacting system. Thus, the 
natural interpretation is that MH-MIT appears first for 
k points at or close to the Fermi surface. Points far away 
from the Fermi surface, on the other hand, like points 4 
and 5 (5, for example, corresponding to F, respectively 
-R), experience only weak renormalizations. The interpre- 
tation is further supported by the influence of non-zero 
t' , which moderately modifies the mass due to changing 
distances of the cluster K points from the Fermi surface 
and a different critical J7c, but otherwise shows a simi- 
lar behavior. We, however, do not observe a significant 
variation of m*j^ across the Fermi surface. The difference 
between points 1 and 2 in Fig. |9] can be explained by 
the large distance of the cluster center from the Fermi 
surface. 

Note that this behavior of the three-dimensional model 
is different from the two-dimensional Hubbard model, 
where the Mott transitio n withi n DCA has been in- 
vestigated in some detail.'^S^IIllI In this case, the so- 
called "momentum selective Mott transition" has been 
observed, where different parts of the Fermi surface un- 
dergo a metal insulator transition for different interac- 
tion strengths. To conclusively exclude the possibility 
of a momentum selective MH-MIT in three dimensions, 
larger clusters with a finer grid of K points on different 
parts of the Fermi surface would be needed. 

In the following we focus on cluster momentum 1 in 
Fig.[9]jb), which is situated directly on the Fermi surface, 
midway between F and M. Since this point exhibits the 
strongest mass enhancement in the insulating phase, it 
is an ideal candidate for studying the MH-MIT. The ef- 
fective mass of this K point is plotted in Fig. 10 as a 
function of U. 

At Uc — 0.766 W, both an insulating and a metallic so- 
lution can be stabilized, depending on the initial Green's 
function used to start the DCA self-consistency. This 
behavior indicates a coexistence region in this regime of 
interaction parameters and tells us that the qualitative 
physical properties of the paramagnetic MH-MIT do not 
change at least qualitatively for a true three-dimensional 
system. The figure also shows the corresponding curves 
for next-nearest-neighbor hopping parameters t' = —0.1 1 
and t' — — 0.2i. Here the coexistence region has vanished 
at the temperature for which the simulations were done, 
while the relatively smooth shape of the curve indicates 
that one is still observing a crossover and not yet a sharp 
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FIG. 10: (Color online) Quasi-particle mass estimate 

at the midpoint between V and M [point 1 in Fig. [9|b)] , for 

T = 0.01 W as function of the interaction parameter U . 



phase transition as in the case of t' 0. This is again 
in accordance with previous DMFT calculations, where 
a reduction of both the critical temperature and the crit- 
ical value of U was observed with increasing t'W^ It is, 
however, different from the interaction transition in the 
two-dimensional Hubbard model, '^'^ where increasing 
leads to an increase of the critical interaction strength. 

It would be highly desirable to perform simulations at 
lower temperatures for non-zero t', but since the com- 
putational effort necessary increases dramatically with 
decreasing temperature, we were not yet able to do these 
simulations for the time being. 

Previous studies of the interaction driven MH-MIT 
at non-zero temperature were largely performed on a 
Bethe lattice in the limit of infinite dimension within 
the DMFT approximation,^* respectively, for a two- 
dimensional Hubbard model using the corr elator projec- 
tion methodic or cluster DMFT.lifiliSlMfUl while the 
general features of the MH-MIT appear to be rather 
insensitive to the actual non-interacting band struc- 
ture, the details like critical values for temperature and 
Coulomb interactions vary strongly with details of the 
model as well as the approximations involved in the com- 
putation and finite size effects. 

To compare values for lattices with different noninter- 
acting density of states yo(w) , BullcP^l suggested using the 
second moment of the DOS 



Won = 4 



Wj2 



dw 



(21) 



-Wl2 



as the characteristic energy scale instead of the band- 
width W. From Eq. (21) one obtains W^s = W for 
the Bethe lattice and WeS ~ 0.816 W for the simple cu- 
bic lattice, and a rather good agreement of critical val- 
ues when relating them to W^s^^^^ Our result Uc = 



0.77(3) W for the coexistence region then translates to 
Uc = 0.94(3) WeS at T = 0.012 Weff. For a conventional 
DMFT calculation, Refs.f^ancPH located the coexistence 
region for this temperature around Uc — 1.18(2) VFoff- 
This indicates that for a true three-dimensional system 
the critical values of the MH-MIT will be renormal- 
ized, and in particular the critical U^ will be shifted 
to lower values. We attribute these renormalizations to 
the short-ranged antiferromagnetic fluctuations present 
in the DCA. They will have the tendency to suppress 
the formation of quasi-particles and will thus cause the 
transition to shift to smaller Coulomb repulsions U. A 
detailed investigation of the location of the transition in 
the thermodynamic limit would require the study of a 
range of cluster sizes and a ca reful finite-size analysis 
along the hues of Refs.ESEIIlIMl. 



IV. ANTIFERROMAGNETIC PHASE 

As previously mentioned the thermodynamically sta- 
ble low-temperature phase of the Hubbard model at half 
filling is, for t' — 0, antiferromagnetic. This phase com- 
pletely covers the MH-MIT. ^■^ Conventionally, the on- 
set of antiferromagnetic long-ranged order is indicated by 
a divergence of the staggered susceptibility upon cool- 
ing from the paramagnetic state at high temperature.!^ 
Complementary to a paramagnetic simulation and analy- 
sis of the susceptibility, the symmetry-broken phase may 
be simulated directly. While this scheme is less accurate 
at determining the location of phase boundaries, it can 
address the properties within the ordered phase and thus 
is relevant for comparison to experiments within that 
phase. Furthermore, it often is desirable to investigate 
the direct change of quantities in the presence of com- 
peting phases or orders. For these reasons we present 
in this paper results obtained in the antiferromagneti- 
cally ordered state. We will show that within the same 
framework, with minor modifications, we can also obtain 
high-quality spectra from QMC data in a phase with non- 
trivial broken symmetry. 

The antiferromagnetic order breaks the translational 
symmetry of the lattice, leading to a doubling of the unit 
cell. This implies that the first Brillouin zone, corre- 
spondingly, halves its size. The resulting magnetic Bril- 



louin zone (MBZ) is shown in Fig. 11 

In order to explicitly break the full translational sym- 
metry along with the SU(2) symmetry, we add a stag- 
gered magnetic field hi = Hq e''^ ''* with Q = (tt, tt, tt) to 
the Hamiltonian of Eq. ([!]) as 



(22) 



where rrii — Ui^^ — rii^^ is the spin polarization at lattice 
site i. In principle this allows the study of properties as 
functions of this staggered field. However, because such 
a field is of little experimental relevance, one is conven- 
tionally only interested in the limit ho ^ 0. If in this 
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FIG. 11: (Color online) Magnetic Brillouin zone of the simple 
cubic lattice and the path along high symmetry points used 
for the presentation of spectra in the antiferromagnetically 
ordered state in Fig. [14] 
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limit a non-zero polarization remains, we have found a 
state with spontaneous symmetry breaking. 

In the actual simulation we add a small field and explic- 
itly break the symmetry (in our case we chose /iq = 0.01) 
in the initialization of our iteration process. The field 
is switched off after the first few iterations and the sys- 
tem is allowed to evolve freely. Eventually, the process 
converges to a solution with either a vanishing staggered 
magnetization Mi{T) oc {rrii) = e'*^ ''' ms(T) = 0, indi- 
cating a parameter regime where the thermodynamically 
stable state is paramagnetic, or ^ and thus an an- 
tiferromagnetically ordered state. 

To be able to include such a field in our simulations, 
we have to ensure that the cluster we use has the proper 
translational symmetry with respect to a double unit cell. 
These clusters are also referred to as bipartite clusters. 
We again employ the systematic classification by Betts'^ 
to find the optimal cluster of this type with 18 sites. Since 
the DCA is formulated in momentum space, the broken 
translational symmetry introduces explicit non-diagonal 
elements in quantities like the Green's function or the 
self-energy. For the following we will adopt the notation 

^ ^ exp [i {Ki R,-K2 Rj)] G„,j {lu^) (23) 

ij 

as extension of Eq. ([s]) for Green's functions with mo- 
menta Ki ^ K2- With this notation, Green's function 
in the antiferromagnetic phase can be represented by the 
2x2 matrix 



K 



(iw„) := 



GaK',K'+Q{^^n) G a K' +Q ,K' +Q(}^n) 



(24) 



where K' is an element of the MBZ. The symmetry re- 
lations G°°^,(iw„) = Gii^.(iw„) = - (Gii^,(iw„))* = 



FIG. 12: (Color online) Staggered magnetization ms{T) as 
function of T for U = 0.671^ as obtained from a DCA calcula- 
tion with Ai' = 18 (circles) and the DMFT (squares). Dashed 
lines: guide to the eye. 



-{GI'k'MY and GiV(i^n) = Gl'^,^ = 
G]^j^,{iuJn) — G°^, (ia;„) hold for Green's function as 
well as for the self-energy. The latter is still defined via 
Dyson's equation 



S<Ti<:'(iw„) = ^c,jc'(i^^«) ^ - Gcrj<:'(iw„)" 



(25) 



which, however, now involves quantities which are 2x2 
matrices. 

This matrix structure makes it necessary to adapt the 
CT-QMC algorithm accordingly. This can most conve- 
niently be done using a spinor representation for the field 
operators, and rewriting the formulas with these new 
composite operators. A detailed account of this proce- 
dure will be given elsewhere. Here we just want to note 
that one can again perform measurements in Matsubara 
space directly and thus obtain an accurate estimate of 
the self-energy, which can then be analytically continued 
as before. 

A first simple test of the method is to calculate the 
staggered moment 



= ^TOi exp(iQ • ri) 



(26) 



and thus locate the antiferromagnetic phase. The re- 
sults of such a calculation for U = 0.67 as function of 
temperature are shown in Fig. [12] for DCA simulations 
(circles) and, for comparison, DMFT calculations using 
Wilson's numerical renormalization group algorithm as 
the impurity solvei'^ (squares). The first thing to note 
is that within DCA the critical temperature is reduced 
by roughly 30% as compared to the DMFT. The val- 
ues of Tn « OmW for DCA and Tat « 0.04214^ for 
DMFT nicely agree with the results obtained by Kent et 
alWi. This is in agreement with the expectation, that 
for a three-dimensional system far enough away from 
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the critical region one should not see dramatic influence 
by the order parameter fluctuations any more. Note, 
however, that to^(O) « 0.39 for both DMFT and DCA 
is reduced as compared to the Hartree approximation, 
where m.sH{Q) k, 0.426. Finally, while the functional 
shape of ms{T) for the DMFT nicely follows the stan- 
dard mean- field behavior m,.(T Tjv) cx y/l - T/Tn 
and, respectively, ms{T — > 0) oc 1 — 2e~^^™/"^, the form 
obtained from DCA is very different, rather exhibiting a 



linear behavior just below T^r and a constant value for 
r< 0.02VF. 

With the ability to perform reliable calculations in the 
symmetry-broken phase, one is of course interested in 
extracting dynamics from the simulations, preferably by 
analytically continuing the single-particle self-energy. To 
this end we again need the high-frequency behavior of the 
self-energy, which can be obtained from a high-frequency 
expansion (see the Appendis]]) as 



icj„ 



(n^) (1 - (n^)) + (m^)^ 



(TOff) {l-2{ns)) 



(mg) (1 - 2{ns)) {ns) (1 - (n^)) + {m-.f 



(27) 



using the staggered spin polarization 



(28) 



The direct analytic continuation of non-diagonal self- 
energies — or Green's functions — is not possible, since the 
non-diagonal spectral function — ^Im E^*]^, (w) has both 
negative and positive values while the standard MEM 
algorithm can only deal with non-negative spectral func- 
tions. In order to solve this problem, we employ the 
linear transformatiort^^ (omitting spin and frequency de- 
pendencies) 



voo 



yll 



(29) 



and determine Im E^^^, {uj) along with the diagonal el- 
ements Im E^*^, (w) 'and Im E^^, (w) using the MEM. 
As in the paramagnetic case, the high-frequency coef- 
ficients of Eq. ( 27 1 are used to normalize the self-energies 



prior to the analytic continuation. Finally, the real parts 
are calculated by a Kramers-Kronig relation analogous 
to Eq. (16). Since the transformation (29) is linear, it 



holds for the analytically continued functions as well and 
can thus be solved for the non-diagonal element E^*]^, (uj). 
An example for a complete self-energy matrix on the real 
frequency axis obtained by this procedure is shown in 
Fig. [TSj Note that the diagonal and off-diagonal elements 
have different symmetry properties. For the particle-hole 
symmetric situation presented here, the former obey the 
relation E;^^,(w + iO+) = E^^,(-w + iO+) and, respec- 
tively, S^^,(a; + iO+) = E;^^,(-a; -|- iO+) following from 
the structure of the Neel state, while the latter are all 
identical but obey T.'^%,{uj + iO+) -E^^,(-w + iO+). 
This last relation in particular implies that the real part 
is an even function of w, and the imaginary part is odd. 

The resulting self-energies for the cluster K' points are 
then interpolated as in the paramagnetic case, and finally 
the spin-averaged spectral function for all momenta k' of 
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FIG. 13: (Color online) Self-energy for K' = V , V = 
0.5 W, il = 0, and T = 0.02 W in the antiferromag- 
netic phase. The real and imaginary parts of the ele- 

" -uj) = E55f'(-t^) and 
^01 I, ,^ _ vio (-, ,\ _ vOi 



ments E°5f,(a;) = 



^00 

are shown. 



the magnetic Brillouin zone follows from 



1. 



AyiiS) = ImTr 



UJ 



(k' 



(30) 



where Tr denotes the trace over the 2x2-matrix. 

Results for single-particle spectra in the antiferromag- 
netically ordered phase and different values of U for 
T = 0.02 W are shown in Figs.[l4l^a)-[l4|c) for paths con- 
necting high symmetry points in the magnetic Brillouin 
zone in Fig. [l3{^b). As expected, the spectral function 
and DOS have a gap around the Fermi energy, i.e., we 
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FIG. 14: (Color online) Spin-averaged single-particle spectra for T = 0.02 W in the antiferromagnetic phase for different 
interaction strengths. Panel (d) shows a spectrum for non-zero next-nearest-neighbor hopping t' = —Q.2t. The left part of 
each figure depicts the local single-particle spectrum for both the majority spins (full red line) and the minority spins (dashed 
blue line). In (a) the edges of the gap are too sharp to resolve properly. To avoid a numerical division by zero, an artificial 
imaginary shift —15 with 8 — 0.03 W was added to the self-energy. The result is a slight broadening of th e ga p edges. The 
interpolation follows the path along the high symmetry points of the reduced Brillouin zone depicted in Fig. 13 b). 



always have an insulating state. Furthermore, it is sym- 
metric with respect to the Fermi energy, reflecting the 
back-folding of the spectrum due to the broken trans- 
lational symmetry. Along a large part of the Brillouin 
zone one has rather flat bands. For weak and moder- 
ate coupling these structures have a rather high spectral 
weight, which results in the formation of characteristic 
van Hove singularities at the gap edges. This is a typ- 
ical weak-coupling result consistent with a conventional 
Hartree approximation. The gap increases for increasing 
[/, while at the same time the weight in the structures at 
the gap edges is redistributed to larger energies, leading 
to a softening of the structures in the DOS. 

Figure [l4[ d) shows an antiferromagnetic spectrum for 
non-zero t' = —Q.2t. The next-nearest-neighbor hop- 
ping breaks the symmetry with respect to the Fermi en- 
ergy analogous to the paramagnetic case. Since magneti- 
cally frustrating interactions cause the antiferromagnetic 
phase to quickly vanish at the present temperature, no 
larger value of t' was simulated. 

V. CONCLUSION 

The dynamical-cluster approximation is a controlled, 
systematically improvable approach to the solution of 



strongly correlated electronic model systems. Due to 
the size and complexity of the cluster impurity prob- 
lem, only quantum Monte Carlo methods are available 
as efficient and unbiased quantum impurity solver algo- 
rithms. While they yield in principle arbitrarily accurate 
imaginary-frequency data, extracting spectral functions 
from QMC data is an ill-posed numerical problem and 
hence remains a difficult task. It requires an analytical 
continuation based on maximum entropy or a similar pro- 
cedure and a thorough error analysis for the quantity to 
be continued including analytically supplemented high- 
frequency information. The previously used Hirsch-Fye 
QMC impurity solver algorithm made the direct contin- 
uation of the irreducible self-energy prohibitively expen- 
sive and made comparatively unreliable root-searching 
techniques necessary. 

The advent of modern continuous-time Monte-Carlo 
algorithms allows a direct simulation of data in frequency 
space and moreover yields high-quality data for the self- 
energy with reliable error estimates, thus allowing a di- 
rect analytical continuation of the self-energy. 

Based on this new route we presented a method to 
extract momentum-resolved dynamical correlation func- 
tions from QMC simulations of strongly correlated elec- 
tron systems. We showed that QMC simulation of clus- 
ters within the dynamical cluster approximation can pro- 
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vide data accurate enough to enable the calculation of 
both momentum- and frequency-resolved single-particle 
spectra. The method can resolve detailed structures in 
the spectral functions, including waterfall-like features. 
We observe that even in three dimensions, momentum re- 
solved self-energies lead to spectra that are qualitatively 
different from dynamical mean-field spectra, and present 
a more reliable starting point for an extrapolation to the 
lattice system. 

In addition we showed that we can access momentum- 
and frequency-resolved spectra in the paramagnetic state 
as well as in more complex ordered phases. As an ex- 
ample we discussed spectral properties of the Hubbard 
model inside the antiferromagnetic phase, also including 
an additional magnetic frustration introduced by a next- 
nearest neighbor hopping. 

We detected the interaction driven metal-insulator 
transition at T = 0.012 Woff and Uc = 0.94(3) PFoff, a 
value which is substantially smaller than the DMFT re- 
sult Uc = 1.18(2) Weff. For non-zero values of t' no MH- 
MIT could be detected for temperatures T > 0.01 W. 
The presence of a smooth cross-over indicates that the 
transition has moved to lower temperatures — an effect, 
that that has been previously studied within the context 
of DMFT.-S2. To clarify this point as well as the t' de- 
pendence of Uc further studies at lower temperatures are 
necessary. Furthermore, in contrast to the standard ex- 
pectation, the MH-MIT appears not to be purely local, 
but rather occurs initially for k points at the Fermi sur- 
face only. This fc-dependent behavior is different from 
the momentum-selective Mott transition found for the 
two-dimensional Hubbard modef^ where the Mott tran- 
sition occurs only on selected parts of the Fermi surface. 
We do not observe such a behavior in our data, but for 
a definite statement on this issue calculations on larger 
clusters would be required. 
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Appendix: High-frequency expansion of the 
self-energies 

We extend the results obtained for the dynamical mean 
field theory^^^^ to the momentum-dependent case en- 
countered in the DCA, both for the paramagnetic and 
the antiferromagnetic phase. 

The antiferromagnetic coarse-grained Green's function 
Go-j<:'(i'^n) ill the presence of a staggered magnetic field 
can be described bj0 



1 

V 



dk 



K/2 



K/2 



K'+k 



(A.l) 



using the matrix notation of Eq. ( 24 1 and S^k' = ^k' ^ 
To gain an expression for the high-frequency coefficients 
of the self-energy, we use the ansatz 

S.K' (ic^n) - Kk' + + O ((ic^„)-') (A.2) 

and expand the coarse-grained Green's function up to 
third order: 



G^j<:'(iw„) 



The result is 



p2 



(-13 



(iw„) (iWn) 



0((i, 



pi 



p2 



p3 



1 
1^ 

ha/2 £,K'+Q 



■'crK' ' 



ha/ 2 S,K'+Q 



VO 
^aK' 



■'aK' 



(A.3) 

(A.4) 
(A.5) 

(A.6) 



where the over-lined quantities are coarse grained over 
the momentum patch centered around , e. g. 



K' 



1 

V 



dk^J,, 



K'+k 



(A.7) 



A direct calculation of the Green's function using 
Heisenberg's equations of motion provides the informa- 
tion necessary for the determination of the unknown coef- 
ficients and Starting again with the Hamil- 
tonian (22), the high frequency coefficients of the single- 



particle Green's function in real space 



C. 



C. 



cl 



Ga^.M = -^ + -^ + -^+0 {M-') (A.8) 



1W„ {iuJnf {iOJnf 



can be obtained^ via 



^aio - 



(A.9) 
(A.IO) 
(A.U) 
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Here [A, B] {{A, B}) denotes the (anti)commutator of the 
operators A and B. A straightforward calculation yields 



^1 



^5. 



(A.12) 

- [K^ + U{nsi)] Sij , (A.13) 



(A.14) 



where S^ij = i + (/i + ^) Sij for nearest neighbors, 



i' + (/i + ^) for next- nearest neighbors, and zero oth- 
erwise. The high frequency coefficients of the coarse- 
grained Green's function in cluster momentum space 
[Eq. (A.3)] are readily calculated by a Fourier trans- 
formation of Eqs. ( A.12 )-( A.14) followed by a coarse- 



graining in fc-space. One obtains 



p2 



'1 



(A.15) 
(A.16) 



tK' 



{U{m^) + K/2) (Ik'+Ik'+q) + U^{m-„) + UK{n^) 



, (A.17) 



where ^k' = £,k' — U/2 and 



^cr) ^ ^ c'^ ij^ai ^ai 



(A.18) 



which is the solution shown in Eq. 27 The non-diagonal 
parts vanish for (too-) — and the expression simplifies 
to the paramagnetic solution of Eqs. ( 11 ) and ( 12 1. 



A comparison with Eqs. A.4 - A. 6| yields 



(A.19) 



1 I {ria) (1 - {ns)) + (mff)2 
(m^) (1-2(^5)) 

(m^)(l-2(ns)) 
{ns) (1 - {ns)) + {msf 



(A.20) 
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